function dy = IB_A(t, y, V_c, V_IB_A,J_a)
%IB_A defines the dynamics of the axon of a deep IB cell.
%   y = [V h_a m_a m_KMa]
%   V_c = [V_sa]
%   V_IB_A = [V_IB_A1 V_IB_A2... V_IB_Ax]

% J_a = -6; (Mark's) % 3; (Aliaa's)
g_NaFa = 100; g_KDRa = 5; g_KMa = 1.5;
g_sa = 0.3; 
g_Leak = 0.25; g_gja = 0.002;
C = 1;

dy = zeros(4,1);

I_Leak = g_Leak*(70+y(1));
I_NaF = g_NaFa*(m_0e(y(1))^3)*y(2)*(y(1)-50);
I_KDR = g_KDRa*(y(3)^4)*(95+y(1));
I_KM = g_KMa*y(4)*(95+y(1));
Is_sa = -1*g_sa*(V_c-y(1));
Ic_gj = 0;
for i = 1:length(V_IB_A)
    Ic_gj = Ic_gj + g_gja*(y(1)-V_IB_A(i));
end

dy(1) = (-1/C)*(J_a+I_Leak+I_NaF+I_KDR+I_KM+Is_sa+Ic_gj);

dy(2) = (1/tau_he(y(1)))*(h_infe(y(1))-y(2));

dy(3) = (1/tau_me(y(1)))*(m_infe(y(1))-y(3));

dy(4) = alpha_KM(y(1))*(1-y(4))-beta_KM(y(1))*y(4);

end  